Plasma metabolomics profiling identifies new predictive biomarkers for disease severity in COVID-19 patients

Recently, numerous studies have reported on different predictive models of disease severity in COVID-19 patients. Herein, we propose a highly predictive model of disease severity by integrating routine laboratory findings and plasma metabolites including cytosine as a potential biomarker of COVID-19 disease severity. One model was developed and internally validated on the basis of ROC-AUC values. The predictive accuracy of the model was 0.996 (95% CI: 0.989 to 1.000) with an optimal cut-off risk score of 3 from among 6 biomarkers including five lab findings (D-dimer, ferritin, neutrophil counts, Hp, and sTfR) and one metabolite (cytosine). The model is of high predictive power, needs a small number of variables that can be acquired at minimal cost and effort, and can be applied independent of non-empirical clinical data. The metabolomics profiling data and the modeling work stemming from it, as presented here, could further explain the cause of COVID-19 disease prognosis and patient management.


Introduction
The SARS-COV-2 pandemic, which has gripped the world over the last three years, has resulted in more than 530 million reported infections and 6.3 million deaths worldwide so far [https://covid19.who.int/]. The pandemic also resulted in unimaginable suffering to individuals, families, communities and countries across the globe. At its peak, the pandemic stressed healthcare systems in different parts of the globe to their limit, disrupted supply chains, destroyed businesses, resulted in massive unemployment and poverty and a never-seen-before upward re-distribution of wealth; which will collectively continue to impact life on earth for possibly generations to come [1][2][3]. Although it is arguable whether the pandemic was July 17, 2020. Patients were diagnosed with COVID-19 using a nasal swab PCR test and later divided into three groups (asymptomatic, mild, and severe) based on their clinical presentation. Each donor gave a 10 mL blood sample, one half of which was collected in a plain tube and the other half in an EDTA vacutainer. A total of 85 samples were collected (30 COVID-19-positive asymptomatic, 10 COVID-19-positive with mild symptoms, and 45 COVID-19-positive with severe symptoms) for the purpose of this study. COVID-19-positive asymptomatic individuals were identified as a result of the national screening campaigns. Symptomatic COVID-19 patients were classified into mild or severe based on guidelines published by Abu Dhabi Department of Health (circular number 33, 19 th April 2020). Patients with mild disease presented with upper respiratory tract infection and symptoms like fever, dry cough, sore throat, runny nose, muscle and joint pains without shortness of breath. Patients with severe disease presented with severe pneumonia and symptoms like fever, cough, dyspnea and fast breathing (>30 per minute), in addition to oxygen saturation <90%. Patient records showed that many of the patients with symptoms were self-medicating with aspirin prior to their hospital visit and that some of the patients with moderate-severe symptoms were placed on dexamethasone and/or heparin subsequent to hospital admission. Immediately upon sample collection, the hospital laboratory staff separated and tested the serum for CRP, D-dimer, ferritin, IL-6 and LDH; a complete blood count was also performed on each sample. Whole blood samples were also aliquoted and frozen at −80˚C for subsequent processing and analysis. The study was jointly approved by the Ministry of Health, Abu Dhabi and Dubai Health Authority (DOH/CVDC/2020/1949) on the understanding that samples will be numbercoded to hide patient's identity, that no personal information will be shared with a third party and that no sample analysis can be performed by entities other than the Research Institute of Medical and Health Sciences (RIMHS), the University of Sharjah (UOS) without prior written approval. No informed consent was required as per the ethical approval decision (DOH/ CVDC/2020/1949); in compliance with the said decision, all samples were fully anonymized before accessing or receiving them.

Serum levels of hepcidin and soluble transferrin receptor (sTfR), sCD163 and haptoglobin (Hp) concentration and phenotype distribution
Upon receipt of frozen samples at RIMHS, UOS laboratories, whole blood samples were thawed and centrifuged; serum was separated and levels of hepcidin (Cat No.733228; MyBiosource, San Diego, California, United States), soluble sTfR (Cat No. 750294; MyBiosource), sCD163 (MBS508555) and Hp (MBS763395) were measured using commercially-available colorimetric assay kits; absorbance was read at 450 nm on a microplate reader. Hp phenotypes were determined by vertical polyacrylamide gel electrophoresis, and the bands were visualized by staining with benzedine solution as previously described [17].

Liquid chromatography tandem mass spectrometry (LC-MS/MS)
Plasma was obtained after the collection of samples into heparinized tubes followed by centrifugation for 5 minutes (3000g). The samples were stored at -80 ºC for long-term storage until further metabolomics analysis. An aliquot of plasma sample was first placed into a microcentrifuge tube where cold methanol was added into the sample at 3:1 v/v (i.e., 30 μL sample, add 90 μL cold methanol) vortex and was then stored at -20ºC for two hrs. Next, the samples were centrifuged at 20,817 x g for 15 min at 4ºC. Then, the supernatant was transferred to a new microcentrifuge tube. Usually, the original sample volume is transferred three times (i.e., for 30 μL sample, add 90 μL cold methanol, then transfer 90 μL supernatant). The sample was dried down using Speed vac at 30-40˚C. The dried sample was then either stored in a -80ºC freezer for further use or dissolved in solvent for LC-MS/MS analysis. Metabolites were analyzed by HPLC-Q-TOF MS/MS using the auto MS/MS positive scan mode as per described in our recent publications [18,19]. Briefly, samples were chromatographically separated using a Hamilton 1 Intensity Solo 2 C18 column (100 mm x 2.1 mm x 1.8 μm) and eluted using 0.1% formic acid in water (A) and 0.1% formic acid in ACN (B) using the following gradient: at a flow rate of 0.250 mL/min 1% B from 0-2 min, then gradient elution to 99% B from 2-17 min, held at 99% B from 17-20 min, then re-equilibrated to 1% B from 20-30 min using a flow rate of 0.350 mL/min. The autosampler temperature was set at 8˚C and the column oven temperature at 35˚C. The ESI source with dry nitrogen gas was 10 L/min, and the drying temperature was equal to 220˚C with nebulizer gas pressure set to 2.2 bar. The capillary voltage of the ESI was 4500 V and the Plate Offset 500 V. MS acquisition scan was set at 20-1300 m/z and the collision energy at 7 eV. Sodium formate was injected as an external calibrant between 0.1 and 0.3 minutes. A total volume of 10 μL sample was injected into the TIMS-TOF MS. Processing analysis was performed using MetaboScape 1 4.0 (Bruker Daltonics). Analyte bucketing and identification were done using the software provided available T-ReX 2D/3D workflow with the following parameters: intensity threshold greater than 1000 counts and peak length equal to 7 spectra or greater. Feature quantitation performed using peak area, for features present in at least 3 (of 12) samples (per cell type) were considered for statistical analysis. Analyte MS 2 spectra were averaged on import and only features eluting between 0.3 and 25 min with mz between 50 and 1000. For metabolite identification, both MS 2 spectra and retention time (RT) were used with the MS/MS spectra as the minimum criterion for a positive hit. For the set of compounds meeting this criterion (either MS/MS alone or MS/MS with RT), annotation using Bruker's implementation of the Human Metabolome Database (HMDB-4.0) was performed; all selected compounds were matched against this library. Where a particular database entry was matched by multiple features, putatively matching features were filtered by considering each of the features against the highest annotation quality score (AQ score) among other putative matches for the same metabolite; i.e. features exhibiting the best fit across the greatest number of factors such as retention time, MS/MS, m/z values, analyte list and spectral library matching were ranked first for the associated identification as per previous publication [18,19]. Pathway enrichment analyses were performed using MetaboAnalyst V5 (https://www.metaboanalyst.ca). Pathway enrichment evaluates overall pathway impact by considering the relative importance of altered metabolite based on their position in the pathway map.

Data analysis
The metabolomics data were first tabulated in Microsoft excel format and then exported to the Statistical Package for Social Sciences (SPSS) software, version 27 [20]. Demographics, clinical and metabolites data were all merged into one SPSS dataset. Descriptive statistics was used to conduct univariate analysis; frequencies and relative frequencies were used to condense categorical data while measures of central tendency were performed for scale data. Normality of scale data was first tested graphically, using Q-Q plots and histograms and then statistically analyzed using the Kolmogorov Smirnov test. Mean and standard deviations (SD) were reported for scale variables showing normally distributed data, whereas median and interquartile range (IQR: Q1-Q3) were used to summarize variables with skewed data. Chi-square test was performed to test for associations between categorical variables where the strength of an association was measured using the odds ratio (OR). To study the relationships between a normally distributed outcome and a categorical dichotomous predictor, the independent t-test was used. If the predictor had more than two groups, one-way ANOVA test was used to compare three or more means. For skewed outcome variables, similar analyses were conducted using the non-parametric tests Mann-Whitney U test and Kruskal-Wallis test, respectively. Spearman correlation coefficient was performed to investigate the correlation between two variables with skewed scale data. P-values less than or equal to 0.05 indicated statistical significance. Bonferroni correction was used to adjust p-values in pairwise comparisons.
The receiver operating characteristic curve (ROC) and the area under the curve (AUC) were performed to identify, from among the clinical and metabolite tests, significant diagnostic biomarkers for predicting the severity of COVID-19 infection. An ROC AUC value above 0.70 indicated moderate to high level of accuracy of prediction. For each test's AUC value, statistical significance was assessed against chance by calculating its 95% Confidence Interval (CI) and associated p-value. For each significant diagnostic test/biomarker showing high/moderate accuracy prediction level (AUC > 0.70), data-driven approach was used to determine the optimal cut-off value, specifically, by maximizing the Youden's index (Sensitivity + Specificity-1) [21]. Next, the sensitivity (SN) and specificity (SP), along with their 95% confidence intervals, were calculated for each diagnostic test.
Optimal cut-off values were used to dichotomize each biomarker into low and high levels. A low level was coded as zero and a high level was coded as 1. After excluding biomarkers that were linearly related, predictors were identified to develop a risk scoring system to define a diagnostic model for COVID-19 severity based on a combination of important biomarkers used as predictors. The risk score was calculated by counting, for each patient, the number of biomarkers that were of high levels. The ROC and the AUC, using Youden's Index, were then used to identify the optimal risk score for predicting the severity of COVID-19. Demographics, clinical and serum metabolite laboratory test results were first compared between the three levels of COVID-19 infection severity (asymptomatic, mild and severe). Preliminary analysis has shown that the asymptomatic and mild groups were comparable on most clinical and metabolite results; no significant differences were observed between the two groups. Accordingly, the two groups (asymptomatic and mild) were clustered into a single group (asymptomatic/mild), which was then compared to the severe COVID-19 group to conduct the analysis reported in this manuscript.

Study population demographics
In this study, we analyzed data pertaining to a total of 85 COVID-19 cases (Fig 1), of whom 35.3% (n = 30) were asymptomatic, 11.8% (n = 10) were mild and 52.9% (n = 45) were severe. Males constituted the majority of patients in this study (84.7%, n = 72) as compared to (15.3%, n = 13) females. Mean age of patients was 42 years (SD = 7.73) with a minimum age value of 27 years and a maximum of 62 years. Age was categorized into two groups where 43.5% (n = 37) were 40 years or younger and 56.5% (n = 48) were older than 40 years (Table 1). Age group and gender distributions were comparable between the two groups (asymptomatic/ mild) vs. severe (Table 1).

Laboratory findings (clinical tests) and severity of COVID-19 disease
In the study sample as a whole, inflammatory markers including the C-reactive protein (CRP) and D-dimer had median values of 18.2 mg/L (Q1-Q3: 5.45-115.49) and 0.60 μg/ml (Q1-Q3: 0.31-2.24) respectively and mean values of lymphocyte and neutrophil counts of 1.45 (SD = 0.72) and 7.56 (SD = 4.13) cells/μL respectively ( Table 1). The majority of the inflammatory markers were found to be significantly higher in the severe COVID-19 group relative to the asymptomatic/mild group. For example, CRP had a median value of 5.90 mg/L in the asymptomatic/mild group and 131.22 mg/L in the severe group (U = 1627.50, p-value<0.001). Similarly, D-dimer had median values of 0.28 μg/ml in the asymptomatic/mild group and 1.27 μg/ml in the severe group (U = 161050, p-value<0.001). While lymphocytes were significantly higher in the asymptomatic/mild group (mean = 1.94 cells/μL) than in the severe group (mean = 1.00 cells/μL; t = 7.880, p-value<0.001), neutrophils were significantly higher in the severe (mean = 9.35 cells/μL) as compared to the asymptomatic/mild group (mean = 5.58 cells/μL; t = -4.773, p-value<0.001) ( Table 1). No significant differences were found between the asymptomatic/mild group vs. the severe COVID-19 group, vis-à-vis median values of serum hepcidin or sCD163. However, Hp and soluble sTfR levels were significantly higher (p-value<0.001) in the severe vs. the asymptomatic/mild group; Hp median values were 138.02 vs. 115.73 ng/ml and sTfR median values were 31.61 vs. 21.46 ng/ml (Table 1).

Metabolomics profiles of COVID-19 patients according to disease severity
To investigate the possibility of identifying serum metabolites that help in studying the prognosis of disease severity, metabolomics profiling of plasma samples from patients with no, mild and severe symptoms was performed. A total of 99 metabolites were measured and compared between the asymptomatic, mild and severe cases. Pairwise comparisons showed comparable results for asymptomatic vs. mild patients, hence the grouping of data obtained from asymptomatic patients and patients with mild symptoms as one "asymptomatic/mild" group; much the same as was done in the previous section. Out of the 99 metabolites (S1 Table), 68 have shown significantly different median values between the asymptomatic/mild group and the severe groups. Of these 68, eight (8) metabolites (K_4_Aminophenol, Acetaminophen glucuronide, Cytosine, Elaidic acide, Glycine, Isobutyric, Paracetamol sulfate and Succinylacetone) were significantly higher in the severe group. Additionally, sixty (60) metabolites showed significantly higher values in the asymptomatic/mild group vs. the severe group (Table 2). Next, we conducted an enrichment analysis of the Biological Process gene ontology terms linked with those metabolites. The enrichment pathway analysis using the "small molecule pathway database (SMPDB)" (available in MetaboAnalyst 5.0 software) revealed that the pathways, that the differentially abundant metabolites were most enriched for included the  citrate cycle, phenylalanine metabolism, phenylalanine, tyrosine, and tryptophan biosynthesis, pantothenate and coenzyme A biosynthesis, tryptophan, glycine, and serine (Fig 2A). Additionally, the same data set produced disease-enriched groups for Hartnup disease, acute seizures, critical illness (major trauma, severe septic shock, or cardiogenic shock), and hyperbaric oxygen exposure when it was searched against the "blood disease signatures database" (available in MetaboAnalyst 5.0 software). As further detailed in the text, the bulk of the diseases or conditions that emerged from this research have symptoms that are consistent with those listed among the most severe COVID-19 cases (Fig 2B).

Predicting COVID-19 disease severity
Predicting the severity of COVID-19 was done at two levels, first by using a single biomarker and then a combination of biomarkers. All clinical and metabolites tests were significantly and strongly associated with the severity of COVID-19. The proportion of patients with severe COVID-19 in the high group of each clinical/metabolites test was significantly higher than that in the low group. The strength of association, as measured by the odds ratio, between disease severity and the clinical and metabolites groups, was lowest for CRP (OR = 4.3, 95% CI: 2.6 to 7.1) and highest for D-dimer (OR = 120.0, 95% CI: 25.1 to 572.9) ( Table 4). A risk scoring system was developed to define a diagnostic model for COVID-19 disease severity based on a combination of important biomarkers used as predictors. All clinical laboratory findings and metabolites that were significantly associated with disease severity were considered as important biomarkers. After excluding biomarkers that were linearly related, and based on the statistical and clinical importance of all identified diagnostic biomarkers, we selected six predictors to conduct the risk-scoring predictive model. This model included the five lab findings (D-dimer, Ferritin, Neutrophils, Hp, and sTfR) and one metabolite (cytosine). A score was calculated for each patient. This score corresponded to the number of biomarkers that were of levels above their respective cut-off values (high group). The accuracy of predicting disease severity in the risk-scoring predictive model was reflected in a highly significant AUC value of 0.996 (95% CI: 0.989 to 1.000) (Fig 3C). The optimal cut-off risk score for the risk-scoring predictive model, as determined by Youden's Index, was 3 with a perfect sensitivity of 100% and a specificity of 92.5% (Table 3). Accordingly, all patients with high levels of at least three of the six predictors would be predicted as developing severe COVID-19 and none of the severe cases would be missed out.   Table).

Discussion
In this study, the concentration of several serum analytes and metabolites was measured in COVID-19 patients with no, mild or severe symptoms. Consistent with numerous previous studies, the serum concentration of analytes that are routinely measured in infected individuals including CRP, ferritin, IL-6, D-dimer, IL-6 and LDH was significantly elevated in patients with severe COVID-19 [22][23][24][25][26]. Also consistent with previous work was the observation that the levels of hepcidin [22,25], sTfR [22,23], and Hp [22,24] were either slightly-moderately elevated or not changed in COVID-19 patients with severe disease. Contrary to the suggestion of Zhou et al [25], our analysis showed that hepcidin is not a true predictor of disease severity in COVID-19 patients. Additionally, while data presented here show that the levels of sCD163 were reduced in severely ill patients, other studies have shown that sCD163 levels increase with disease severity [27]. This discrepancy could be a reflection of variations in methodology, sample collection timing, and/or differences in macrophage activity [28]. Irrespective of these discrepancies, variations in sCD163 concentration seem to have little, if any, impact on COVID-19 disease severity. Our data also showed that Hp phenotype distribution was similar in severe vs. asymptomatic/mild groups, which is in agreement with previous work which has suggested that Hp phenotype has no bearing on COVID-19 disease severity [29]. With regard to the metabolomics profiling of COVID-19 patients and as noted earlier, 60 metabolites decreased in the severe cases relative to asymptomatic/mild patients' group. The list of identified metabolites included several amino acids, vitamins and few fatty acids (Table 1). These are regarded as the fundamental elements that support the rise in cellular demands during illness. However, through catabolism pathways, they are also involved in innate and adaptive immune responses to infection [30,31]. Therefore, the decrease in some of the reported metabolites is consistent with previous studies, that reported lower levels of amino acids in hospitalized COVID-19 patients compared to asymptomatic ones [32,33]. The outcomes do in fact support the previously noted negative correlation between amino acids and immune responsiveness and hyper-inflammation indicators [34], which is characteristic of severe COVID 19. For example, the levels of Kynurenic acid in severe cases was found to be lower than that in asymptomatic/mild cases. Previous studies have suggested that the Tryptophan catabolite/ Kynurenine pathway may play a significant role in COVID-19 and critical COVID-19 [35]. Moreover, it appears that the increased level of kynurenine and the ratio of kynurenine to tryptophan is strongly correlated with the severity of COVID-19 patients [32,36,37]. Interestingly, the ratio of Kynurenic acid/Kynurenine did not significantly differ between COVID-19 patients compared with non-COVID-19 controls, indicating no significant changes in Kynurenic acid activity, according to a systematic review [35]. This is indeed consistent with our finding that in patients with severe COVID-19, tryptophan and Kynurenic acid levels were significantly lower than in the counter group (Table 3). Over the last three years, much time and effort has been spent on identifying serum biomarkers that could predict disease severity in COVID-19 patients with high accuracy. Elevated levels of serum biomarkers such as ferritin, IL-6, D-dimer and lactate dehydrogenase among others were reported to be valuable predictors of disease severity and death [22][23][24][25][26]. However, not all COVID-19 patients showing elevated levels in one or more of these biomarkers ended up with severe disease and death [38]. In other words, relying on one or more serum analytes tends to yield low prediction accuracy as evidenced by the fact that such approaches could only account for only a significant percentage of cases. In this context, numerous predictive models that relied on overlapping sets of variables drawn from COVID-19 patients' demographic data, clinical signs and symptoms, chest X-ray imaging and co-morbidities were proposed (reviewed in [39]. However, many of these severity predictive models suffer from a high degree of subjectivity and high likelihood of bias [39].  [40]. Besides the fact that this model could be skewed by the demographics components making it more population-specific than desired, achieving 93% accuracy by relying on 6 variables is cumbersome and difficult to apply in many poor countries and rural settings. Another AI-based model was developed by relying on laboratory findings including LDH, IL-6, D-dimer, fibrinogen, glucose, monokine induced by gamma interferon (MIG) and macrophage derived cytokine (MDC) levels in 60 COVID-19 Russian patients [41]. The model described by the study relied on eight parameters (creatinine, glucose, monocyte number, fibrinogen, MDC, MIG and CRP) to yield a predictive power of 83−87% [41]. In other words, this laboratory findings-based model failed to account for 13-17% of patients at risk of severe disease.
With the availability of six metabolite predictive biomarkers at our disposal, we sought to develop a high accuracy prediction model based on disparate data-integrating approaches; namely, patient laboratory findings and plasma metabolomics profiles. The statistical model developed and tested was based on ROC-AUC values. Although some metabolomes including acetaminophen, acetaminophen glucuronide and paracetamol sulfate significantly predicted COVID-19 severity, it was unlikely that these metabolomes were involved in the pathophysiology of COVID-19. Severe cases of the disease were more likely to receive more paracetamol than asymptomatic/mild cases. Therefore, these metabolomes were excluded from the predictive models. Accordingly, the predictive model was developed, and its prediction accuracy was internally validated using six biomarkers that were not linearly-related; namely, D-dimer, ferritin, neutrophil counts, Hp, sTfR and cytosine. Prediction accuracy of disease severity using this model was 0.996 (95% CI: 0.989 to 1.000) with optimal cut-off risk score of three biomarkers. In other words, out of 100 patients with severe COVID-19 showing significant elevation in at least three of the six metabolites would predict disease severity in all 100 patients. The model has the advantage of yielding high predictive power with small set of variables (three laboratory findings) that can be easily and quickly acquired at minimal cost. Moreover, the predictive model can be dynamically applied independent of non-empirical clinical data (co-morbidities, signs and symptoms and loss of taste or smell among others) and can be dynamically applied as the disease progresses making timely and proper clinical interventions possible. That said, the utility of the model remains limited by the fact that the study was conducted retrospectively on a small number of samples. Another limitation in our study is that, with the number of COVID-19 cases gradually dwindling to almost zero in the UAE as in most parts of the world, we were not able to compile a new independent dataset with the same set of predictors as means of validating our prediction model. Future studies are recommended to test the validity of the suggested model on multiple datasets to ensure its generalizability.

Conclusion
By integrating laboratory findings and metabolomic profiling data, a model to predict disease severity in COVID-19 patients was generated. The accuracy of the model was high (>98%), and it has the advantage of requiring three biomarkers to yield high sensitivity and specificity in predicting disease severity. The suggested model may prove useful in better managing COVID-19 patients at high risk of severe disease. Lastly, the fact that the model included cytosine as a biomarker and that cytosine is not usually included in routine laboratory testing for COVID-19 patients, merit further work on developing reliable and highly sensitive, yet quick and easy to perform, assays to measure serum cytosine concentration.
Supporting information S1